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We investigate the characteristics of two dimensional melting in simple atomic systems via 
isobaric-isothermal (NPT) and isochoric-isothermal (NVT) molecular dynamics simulations with 
special focus on the effect of the range of the potential on the melting. We find that the system 
with interatomic potential of longer range clearly exhibits a region (in the PT plane) of (thermody- 
namically) stable hexatic phase. On the other hand, the one with shorter range potential exhibits a 
first-order melting transition both in NPT and NVT ensembles. Melting of the system with inter- 
mediate range potential shows a hexatic-like feature near the melting transition in NVT ensemble, 
but it undergoes an unstable hexatic-like phase during melting process in NPT ensemble, which 
implies existence of a weakly first order transition. The overall features represent a crossover from a 
continuous melting transition in the cases of longer-ranged potential to a discontinuous (first order) 
one in the systems with shorter and intermediate ranged potential. We also calculate the Binder 
cumulants as well as the susceptibility of the bond-orientational order parameter. 

PACS numbers: 64.70.Dv, 02.70.Ns, 05.70.Fh, 61.20.Ja 



I. INTRODUCTION 

For decades, two dimensional melting|l| has been 
an important subject of research in condensed matter 
physics both theoretically and experimentally. One of 
the most important theoretical frameworks was given by 
Halperin and Nelson 0, and Young Q who proposed 
(building upon the work by Kosterlitz and ThoulessQ) 
the so called KTHNY theory that the two dimensional 
melting can occur in two stages of continuous defect- 
mediated transitions with the intermediate hexatic phase 
characterized by quasi-long-range orientational order and 
short range translational order [y]. 

One of the most important questions in the two di- 
mensional melting, which has not been satisfactorily 
answered yet, is probably the question of how to de- 
termine the form of the inter-particle potential that is 
most favorable for the existence of the hexatic phase. 
Even though several experimental studies support the 
existence of hexatic phases [7H17j|. computational stud- 
ies of two-dimensional melting of hard-core potential 
systems [l8l. [27M30I ] including hard discs or Lennard- Jones 
(LJ) potentials tend to favor first order transition scenar- 
ios (though some conflicting results also exist) @. Chui, 
et al jl9H22j advanced the possibility of first order melt- 
ing transition through grain boundary formations when 
the defect core energy becomes low enough. Kleinert 
and Janke argued that the nature of the two-dimensional 
melting can change from continuous to first order tran- 
sition as the magnitude of the so-called angular stiffness 
of the local rotation field 23H25| is decreased. From this 
argument, they contended that the melting of LJ systems 
would occur via first order, while, it would be continu- 
ous in the case of Wigner crystals where the particles are 
interacting via long-range Coulomb potentials. 

Recently we investigated on the criterion for the exis- 



tence of the hexatic phase by tuning the form of the in- 
terparticle potential (which is the Morse potential) that 
could change the size of the range of the dominant in- 
terparticle interaction (26|. The Morse potential can be 
written as 



V M (r) = 



eo 
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eo 
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where, r is the distance between particles, a the distance 
from the origin to the minimum of the potential, and eo 
is the strength of the interaction. After setting the eo = 1 
and <7 = 1, we can vary the value of the single parameter 
a to tune the softness and the range of the potential. 
Smaller value of a corresponds to a softer potential with 
longer range of the attractive part of the potential. On 
the other hand, as the value of a increases, the potential 
gets stiffer and shorter-ranged (Fig. Q]). 

In our previous work, we investigated the trend of hex- 
atic phase formations as the range of the potential is var- 
ied. Detailed simulation results were presented especially 
in the regime of softer potential with a = 3.0, where the 
melting exhibits a stable region of hexatic phase on th 
PT plane. 

Here in this work, we investigate how the character- 
istics of melting evolve as the value of a is varied from 
a = 3.5 to larger values of a — 6 and a = 12. See Fig. [TJ 
for the shape of the potentials for different values of a. 
We observed that different types of melting are exhib- 
ited for an atomic system described by the Morse poten- 
tial. A system with a = 3.5 exhibits a melting transition 
which is compatible with the KTHNY theory, showing 
thermodynamically stable hexatic phase as well as two 
stage melting. For simulations with a = 6, it is found 
that, although some hexatic-like features are revealed in 
NVT ensemble simulations, it exhibits a weakly first or- 
der transition in NPT ensemble. Systems with a short 
range potential with a = 12 show a strong first order 
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transition compared to the case of a = 6, clearly show- 
ing coexistence of solid and liquid. These results appear 
to be consistent with the arguments of Kleinert[25( in 
that the systems with short-ranged potential are corre- 
lated with smaller angular stiffness and first order melt- 
ing transition. 



II. SIMULATION METHODS AND RESULTS 

In this work, we performed NPT MD simulations us- 
ing the modified Parrinello- Rahman (PR) method[3lll32j| 
combined with Nose- Hoover (NH) thermostat (33|. As for 
the mass of the particles, for convenience, we put m = 1 
which implies that the time unit to = \J ma 2 /eo also 
becomes unity when we set a = 1 and eo = 1. The 
equations of motion were integrated via the Nordsieck- 
Gear 5th-order predictor-corrector method with the in- 
tegration time step of At = 0.002. This guarantees the 
conservation of the total Hamiltonian without noticeable 
drift. In the simulations we used two empirical param- 
eters, the barostat mass M v = 1 and the thermostat 
mass M s = 1. Test simulations with several other val- 
ues (M v = 0.1, 1,10, M s = 0.1,1,10) of the parameters 
were also performed with almost the same results. The 
number of particles employed ranges from N — 400 up 
to N = 10000. 

In order to investigate the characteristics of the melt- 
ing transition, we obtained the isothermal equation of 
state on the plane of pressure vs. density. This was ob- 
tained by the NH-MD simulations by decoupling the PR 
(isobaric) part from NH-PR MD equations of motions by 
taking M v — oo which reduces the system to NVT con- 
dition. The pressure was evaluated by means of the virial 
expression for the range of the densities corresponding to 
the region of transition from liquid to solid. For each 
density, 10 6 ~ 3 x 10 6 steps of integration were carried 
out for equilibration beginning with a configuration of 
triangular lattice and, after equilibration, 10 7 steps of 
integration were performed for thermodynamic calcula- 
tions. 

In this case of NVT ensemble we have to fix the shape 
of the box. In order to reduce the finite size effect, we 
used a rhombic box (with the smaller side angle of 60 de- 
grees) for the shape of the system with periodic boundary 
conditions. However, independent results of ours from 
square box showed no significant difference (from those 
of rhombic box) with respect to the quantities of our in- 
terest. 

Important criterion for the existence of hexatic phase 
(and hence continuous melting transition) would be that 
the isothermal equation of state for pressure vs. den- 
sity exhibit a monotonic behavior together with a non- 
monotonic region of dip in the slope dPj dp. On the other 
hand, a first order melting transition would be associ- 
ated with the existence of a van der Waals type loop in 
pressure vs. density curve with unstable and metastablc 
region in NVT ensemble simulations. 



In order to investigate the nature of the possible hex- 
atic phases, one has to compute the bond-orientational 
order parameter. The local bond-orientational order pa- 
rameter ip6( r ) at position r is defined as 

Mr) = ±J2e™^\ (2) 

1 3 

Here, the sum on particle j is over the Ni neighbors of 
the particle i (corresponding to r at the center with fty 
being the angle between the particles i and j with re- 
spect to a fixed reference axis. We regarded the particles 
within a cutoff radius as the neighbors, where the cutoff 
radius is chosen as the first minimum of the pair correla- 
tion function of the system. This method is found to be 
efficient and reliable for large scale simulations [34j]. 

Then the global bond-orientational order parameter is 
defined as 




r 



where N denotes the total number of particles in the 
system. In order to distinguish the bond-orientational 
order of the different thermodynamic phases, we com- 
pute the spatial correlation function Ge(r) of the bond- 
orientational order parameter, defined as [2j 

G e (r)=<MrW 6 (0)>, (4) 

In the hexatic phase, according to KTHNY theory, the 
bond-orientational correlation function is expected to ex- 
hibit an algebraic decay i.e., Gq(t) ~ r~ Ve with the decay 
exponent j]q < 1/4, where t]q = 1/4 corresponds to the 
limit of the power-law decay behavior in the KTHNY 
theory. 

In order to further understand the nature of the hexatic 
phase we obtained the histogram distribution (35j of the 
density order parameter for different system sizes for the 
values of pressure and temperature where a melting tran- 
sition is expected to occur (from other measurements). 
We expect that histograms with single peaks would im- 
ply that there exists continuous melting transition. On 
the other hand, existence of double peaks with increasing 
peak heights (as the system size increases) would imply 
a first order transition. 

We also investigate the behavior of the linear suscepti- 
bility for the global bond-orientational order parameter 
near the melting transition in order to check the con- 
sistency with the result from the isothermal equation of 
states. Specifically, we obtain the size dependence of the 
susceptibility by calculating the fluctuation of the bond- 
orientational order parameter for sub-blocks of the sys- 
tem with linear sizes L, which is defined as [36] 

XL = L d (<**) £ - (* 6 } 2 L ) . (5) 

where d = 2 is the spatial dimension. In the computation, 
the system is sub-divided into sub-blocks of linear sizes 
with L = N /Mi, where ranges from 10 to 20. 
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Also, it is useful to calculate the Binder cumulant for 
the global bond-orientationl order parameter for subsys- 
tem (linear) size L is defined by 



where the subscript L denotes that the quantities are 
calculates for subsystem sizes of linear size L. 



A. The case of a soft and long-ranged Morse 
potential: a — 3.5 

Here, we first deal with the case of a moderately soft 
(and longer ranged) potential with a = 3.5. Figure [21 
shows the isothermal equation of state (at T — 0.7) on 
the plane of pressure vs. density. This was obtained by 
the NH-MD simulations by decoupling the PR (isobaric) 
part from NH-PR MD equations of motions by taking 
M v — oo which reduces to NVT condition. We here de- 
fine the density p as p = No 1 jV where N is the total 
number of particles and V the total volume (area in two 
dimensions) of the system. Densities were chosen from 
the range p = 1.56 ~ 1.62, with the density increment 
of Ap = 0.005 and the pressure was evaluated by means 
of the virial expression (with ks = 1). This range of 
the density corresponds to the region of transition from 
liquid to solid. For each density, 10 6 ~ 3 x 10 6 steps 
of integration were carried out for equilibration begin- 
ning with a configuration of triangular lattice and, after 
equilibration, 10 7 steps of integration were performed for 
thermodynamic calculations. 

The number of particles employed was N — 3600. In 
order to reduce the finite size boundary effect, we used a 
rhombic box (with the smaller side angle of 60 degrees) 
for the shape of the system with periodic boundary condi- 
tions. However, independent results of ours from square 
box showed no significant difference (from those of rhom- 
bic box) with respect to the quantities of our interest. 

The isothermal curve increases monotonically near the 
transition region satisfying the condition of mechanical 
stability (unlike the discontinuity of density in a first or- 
der transition) that the isothermal compressibility should 
be positive Kt = (l/p)(dp/dP)T > 0. We may identify 
the boundary of stable hexatic phase as the values of the 
density where an abrupt change in the isothermal com- 
pressibility occurs. In this way, we estimate the density 
of solid-hexatic transition as p s -h — 1-6. 

Although the change in isothermal compressibility is 
less conspicuous near the hexatic-liquid boundary, we see 
that, near the density 1.58 < p < 1.585, there exists a 
crossover in the slope of the isothermal compressibility. 
Below, we give an estimation of the density of hexatic- 
liquid transition by applying a theoretical expectation 
from KTHNY theory on the decay exponent of the spa- 
tial correlation of the orientational order parameter (see 
below) . 



The fact that the pressure within the hexatic phase is 
monotonically increasing as the density increases (with 
the resulting isothermal compressibility kept positive) 
appears to be a very compelling evidence for a stable 
hexatic phase in thermal equilibrium. 

In order to distinguish the orientational order of the 
phases, we have computed the bond-orientational corre- 
lation function G&(r) defined above. 

Figure. [3] is the bond-orientational correlation function 
for the range of the density (1.56 < p < 1.615). We find 
that, for the density range of 1.585 < p < 1.595, the aver- 
aged correlation functions exhibit algebraic decays with 
the decay exponent rj < 1/4 while, at lower densities, 
they exhibit decays with faster than power law behav- 
ior in the long distance limit. At p — 1.585, the ori- 
entational correlation exhibit a slope of approximately 
r\ — 0.25. Here, the crossover from a power law decay 
to exponential takes places with exponent approximately 
equal to 1/4, and also this value of the crossover density 
agrees almost precisely with the density region exhibiting 
an abrupt change of the slope i.e., of the compressiblity 
in the equation of state (Figure [2]). 

Now, the fourth order Binder cumulant for the global 
bond-orientationl order parameter for sub-block systems 
of (linear) size L is shown in Fig. 2) We can see that the 
density at the crossing point is around p ~ 1.585. This 
is compatible with the boundary density (1.585) between 
the liquid and the hexatic phase which was shown above 
in isothermal curve for NVT ensemble, and also compat- 
ible with the the boundary density between liquid and 
hexatic-like phase in terms of the power law decay of the 
bond-orientational order near p = 1.585 

It may be expected theoretically that the binder cumu- 
lants of local orientational order in hexatic phase collapse 
to a line because of the critical charateristic of the phase. 
However, we may not consider non-collapse of the Binder 
cumulants to a line as an evidence for non-existence of 
hexatic phases. This is because, even for the case of XY 
model, complete collapse was not found in the region 
where orientational order decay algebraically, but rather 
it exhibits a crossing point at the transition temperature 

ma 

Now, we turn to the linear susceptibility for the global 
bond-orientational order parameters near the melting 
transition. Shown in Fig. [5] is the sub-block susceptibility 
obtained from the fluctuation of the bond-orientational 
order parameter for sub-blocks of the system with linear 
sizes L with L = N/M^ where Mb ranges from 11 to 20. 
We see that the suceptibility shows a broad peak region 
near the density 1.580 < p < 1.585 which borders the 
liquid-hexatic phase boundary region. We also see that 
the suceptibility exhibits broader shape (showing weaker 
dependence on the density) in the liquid region compared 
with other cases that will be shown below. 

Figure [6] shows a snapshot of the particle configuration 
for density p = 1.585 within the hexatic phase region but 
close to the transition (to liquid phase) at the tempera- 
ture T = 0.7, which shows free dislocations (i.e., bounded 
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pair disclinations) . This shows rather clearly the funda- 
mental role of defects leading to the power law decay of 
orientational correlations. 

In order to further understand the nature of the hex- 
atic phase we obtained the histogram distribution 35] 
of the density order parameter for five different system 
sizes (N = 900, 1600,3600, 10000) under constant exter- 
nal pressure and temperature of T = 0.7, and P — 13.5, 
where a hexatic phase is expected to occur from our mea- 
surement of the orientational correlations. In Fig. [7] we 
see that all the histograms exhibit single peaks, which 
implies that there exist a unique phase with minimum 
free energy. It is also observed that, as the number of 
particles increases, the peak height becomes higher and, 
at the same time, the width of the peak decreases. Also 
the position of the peak tends to shift to the lower density 
(within the hexatic regime) probably due to the develop- 
ment of long range fluctuation. This indicates that this 
region corresponds to a single phase region (unlike solid- 
liquid mixture) consistent with the absence of van der 
Waals loop in the pressure vs. density curve. 

All of these observations lead us to the conclusion that 
there exist a thermodynamcally stable hexatic phase con- 
sistent with the KTHNY melting scenario for the case of 
a = 3.5. 



B. The case of an intermediate-ranged potential: 

a — 6 

Next, we examine an intermediate-range potential with 
a = 6, which corresponds approximately to the famous 
LJ potential[38j]. In Fig. [51 the equation of state ex- 
hibits a weak van der Waals-like loop in the pressure vs. 
density, which indicates a first order transition. The un- 
stable region ranges from p ~ 1.04 to p ~ 1.065. This 
is confirmed more rigorously by the histogram distribu- 
tions of the density in NPT ensemble simulations for 
different system sizes. In Fig. |H1 the histogram distribu- 
tions of the density for systems with N = 900, 3600, and 
10000 are shown from NVT ensemble for T = 0.57 and 
P = 1.85. For the cases of N = 900 and 3600, we observe 
transitions between two peaks (through the valley of fi- 
nite height between the peaks) with the resulting double 
peaks in the histograms. For the N = 10000 systems, 
however, we can no longer observe crossing between the 
the coexisting phases. Instead, we could observe two dif- 
ferent (separate) histograms that are determined by the 
initial states, depending whether the initial state is in the 
ordered solid phase or in the disordered liquid phase. Ev- 
idently, the free energy barrier increases with increasing 
system size, which indicates clearly that the transition is 
of first order (3f|. Nevertheless, the system configurations 
in the coexistence region resembles those of the hexatic 
phases (Fig. [T0j) , showing algebraically decaying orienta- 
tional order (Fig fTTT) . We see that the boundary between 
liquid and hexatic-like phase in terms of the orientational 
order is located around p = 1.05 ~ 1.055. 



Furthermore, we also observe a hexatic-like feature in 
the NPT ensemble, where the system goes through tem- 
porarily a hexatic-like phase before transiting into the 
other phase. This kind of characteristics in NPT ensem- 
ble seems to have been already reported as 'metastable 
hexatic phase' in LJ system by Chen et al [39j - l4l1 | . Al- 
though they argued that large system size (N ~ 40000) 
is necessary to observe this kind of features, we could ob- 
serve such metastable hexatic phases even for systems 
with smaller sizes of N = 1600. We think that this 
hexatic-like feature in a first order transition can be at- 
tributed to weakly first order nature of the transition. As 
shown below, the relative free energy barrier in this case 
is considerably lower than that in the case of shorter- 
ranged potential of a = 12, from which we presume that 
the metastable or unstable hexatic-like phase comes from 
defect proliferation by thermal fuctuations, but not from 
some mechanism leading to a true phase transition. 

Now, the fourth order Binder cumulant for the global 
bond-orientationl order parameter for sub-block systems 
of (linear) size L is shown in Fig. [T2j We can see that 
the density at the crossing point is located near p ~ 1.06. 
This corresponds to a point in the middle of the unsta- 
ble part of the van der Waals curve in the isothermal 
equation of state. 

Next, we deal with the linear susceptibility for the 
bond-oricntational order parameters. Shown in Fig. ll3l is 
the sub-block susceptibility obtained from the fluctuation 
of the bond-orientational order parameter for sub-blocks 
of the system with linear sizes L with L = N/Mf, where 
Mb ranges from 10 to 20. We see that the suceptibil- 
ity shows a sharper peak (as compared with the case of 
a = 3.5) at the density p ~ 1.045 which is a little bit 
below the density {p = 1.05 ~ 1.055) where the orien- 
tational correlation exhibits a spatial decay exponent of 
0.25. This might be interpreted as a small evidence that 
the nature of the melting transition of this system is in- 
consistent with the expectation of the KTHNY theory. 



C. The case of a short-ranged potential: a = 12 

Finally, we investigate the case of a short-ranged po- 
tential with a = 12. Figure. [T3] shows the equation of 
state in the density region of 0.85 < p < 1.075 obtained 
from NVT ensemble simulations at T — 0.57. We can see 
that the equation of state exhibits a van der Waals-like 
region in the density with the unstable region ranging 
from p ~ 0.96 to p ~ 1.04 clearly indicating a first order 
melting transition. 

The fourth order Binder cumulant for the global bond- 
orientationl order parameter for sub-block system size L 
is shown in Fig. [15] We can see that the density at the 
crossing point is located near p ~ 1.045 which is located 
near the lower density limit of the metastable solid (spin- 
odal) as shown in the NVT isothermal equation of state. 

First order nature of the melting transition is con- 
firmed further with the double peak nature of the his- 
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togram distributions of the density from NPT ensemble 
simulations for system sizes of N = 400, 900, 3600, 10000 
as shown in Fig. [TH In this case, the free energy barrier 
is presumably much higher than the case of a = 6, and 
none of the systems with a = 12 exhibit any tunneling 
transitions between liquid-like states to solid-like states 
during 10 8 MD steps of simulations with AT = 0.002. 
Therefore, double peak histogram distributions for each 
of the system sizes are actually obtained by combining 
two separate histograms, one with ordered initial states 
(higher density) and the other with disordered initial 
states (lower density), respectively. 

Also, typical system configuration for a — 12 is shown 
in Fig. [17|for p = 1.0 and T = 0.57 corresponding to the 
coexisting region, where we can see that hexatic-like fea- 
ture disappears, and that liquid phase region consisting 
of defect clusters coexists with solid region. Also in NPT 
ensemble simulations of the melting process, we can no 
longer observe metastable or unstable hexatic phase, but 
observe a discontinuous abrupt change in density. From 
these observations we thus conclude that first order na- 
ture of the melting gets stronger as the potential range 
decreases. 

Next, we look into the linear susceptibility for the 
bond-orientational order parameters. Shown in Fig.[T51is 
the sub-block susceptibility obtained from the fluctuation 
of the bond-orientational order parameter for sub-blocks 
of the system with linear sizes L with L = N/Mj, where 
Mb ranges from 12 to 20. We see that the suceptibility 
shows a peak at the density p ~ 0.96 which is located near 
the limit of the metastable liquid (spinodal) as shown in 



the curve of the isothermal equation of state from NVT 
ensemble. 



III. SUMMARY AND DISCUSSIONS 

In conclusion, we have reported on some details of two 
dimensional melting in systems of particles interacting 
via Morse potential when the range of the potential is var- 
ied. We showed that the melting of system with longer- 
ranged potential (a = 3.5) clearly exhibits features of 
melting consistent with KTHNY theory exhibiting sta- 
ble hexatic phases. As the range of the potential de- 
creases, however, we observe a crossover in the transition 
nature to a first order transtkm. In the case of a = 6 
where the range of the potential is intermediate, the sys- 
tem exhibits a weakly first order melting transition with 
some unstable hexatic-like phase during melting process 
in NPT simulations. In the case of a = 12 where the 
range of the potential is shorter, we observe a stronger 
first order melting. It appears that the crossover from 
continuous to first order melting transition in this sys- 
tem is related to the decrease of the so-called angular 
stiffness of the rotation field [23T4251 ] . It would be inter- 
esting to carry out a detailed calculation of the angular 
stiffness in our model system as the value of a is varied. 
It would be also interesting to find a possible connection 
to the change from continuous to first order transition in 
two dimensional XY model when the shape of the XY 
potential gets sharpened p^ - (43 | . 
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FIG. 1: Shape of the Morse potential for a = 3.5, 6 and 12. 
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FIG. 2: Isothermal equation of state (pressure vs. density) 
for a — 3.5 at the temperature T = 0.7 obtained from NVT 
ensemble with N — 3600 particles (see the text for details) . 
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FIG. 3: Spatial correlation functions (for a — 3.5) of the 
bond-orientational order parameter for different densities at 
T = 0.7. The dashed line indicates a power law decay with 
the decay exponent of 1/4. 
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FIG. 4: Fourth order Binder cumulant of the bond- 
orientational order parameter vs. density obtained by sub- 
block method at T = 0.7 (a = 3.5) obtained from NVT 
ensemble. The total number of particles is N = 10000 
(100 x 100 and the sub-blocks have dimensions of L x L = 
(N/Mb) x (N/M b ) with M b ranging from 11 to 20. 
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FIG. 5: Orientational susceptibility vs. density, obtained 
from the fluctuation of the bond-orientational order parame- 
ter using sub-block method with the same system as in Fig. [4] 
at T = 0.7 (a = 3.5). 
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FIG. 6: (Color online) A snapshot of the configuration of par- 
ticles at density p = 1.585 for a = 3.5 (T = 0.7) represented 
with Delaunay triangulation. Blue open circles and red open 
squares denote the defect sites of particles with five nearest 
neighbors and seven nearest neighbors, respectively. 
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FIG. 7: Histogram distributions (for a = 3.5) of the particle 
density from NPT ensemble simulations at P = 13.5, T = 0.7 
for system sizes N = 900, 1600, 3600, and 10000, respectively. 
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FIG. 8: Isothermal equation of state (pressure vs. density) 
obtained from NVT ensemble with N = 3600 (for a = 6) at 
temperature T = 0.57 which shows a weak van der Waals-like 
loop. 
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FIG. 9: Histogram distributions of the particle density (for 
a = 6) from isobaric-isothermal ensemble simulations at 
T = 0.57, P = 1.85 for system sizes N = 900, 3600, and 
10000. Note that we obtain two separate histograms for the 
system size N — 10000 with ordered initial states (higher 
density) and with disordered initial states (lower density) , re- 
spectively. This is because the system of N = 10000 does 
not exhibit transitions between liquid-like states to solid-like 
states during 10 8 MD steps with AT = 0.002. 
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FIG. 10: (Color online) A snapshot of the configuration of 
particles at density p = 1.06 for a — 6 (T = 0.57) represented 
with Delaunay triangulation. The defects are indicated with 
the same method as in Fig. [6] 
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FIG. 11: Spatial correlation functions of the bond- 
orientational order parameter for different densities (for a = 
6). The dashed line indicates a power law decay with the 
decay exponent of 1/4. 
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FIG. 12: Fourth order Binder cumulant of the bond- 
orientational order parameter vs. density obtained by sub- 
block method at T = 0.57 (for a — 6) The total number of 
particles is N — 10000 (100 x 100 and the sub-blocks have 
dimensions of L x L — (N/Mb) x (N/M b ) with M b ranging 
from 10 to 20. 




1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 



Density 

FIG. 13: Orientational susceptibility obtained from the fluc- 
tuation of the bond-orientational order parameter using sub- 
block method with the same system as in Fig.[T2]at T = 0.57 
(for a = 6). 
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FIG. 14: Isothermal equation of state (pressure vs. density) 
obtained from NVT ensemble with N = 3600 (for a = 12) at 
temperature T = 0.57 exhibiting a clear van der Waals-like 
loop. 
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FIG. 15: Fourth order Binder cumulant of the bond- 
orientational order parameter vs. density obtained by sub- 
block method at T = 0.57 (for a = 12) The total number of 
particles is N = 10000 (100 x 100 and the sub-blocks have 
dimensions of L x L = (N/M b ) x (N/M b ) with M b ranging 
from 12 to 20. 
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FIG. 16: Histogram distributions of the particle density (for 
a — 12) from isobaric-isothermal ensemble simulations at 
T = 0.57, P = 1.825 for system sizes N = 400, 900, 3600, 
and 10000. Note that none of the systems exhibit any tunnel- 
ing transitions between liquid-like states to solid-like states 
during 10 s MD steps of simulations with AT = 0.002. There- 
fore, we obtained two separate histograms for each of the 
system sizes with ordered initial states (higher density) and 
with disordered initial states (lower density), respectively. 




FIG. 17: (Color online) A snapshot of the system configura- 
tion at density p = 1.0 for a = 12 (T = 0.57) represented 
with Delaunay triangulation. The defects are indicated with 
the same method as in Fig. [6] We can clearly see coexistence 
of liquid-like region and solid-like region. 
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FIG. 18: Orientational susceptibility obtained from the fluc- 
tuation of the bond-orientational order parameter using sub- 
block method at T = 0.57 (for a — 12). The total number 
of particles is N = 10000 (100 x 100 and the sub-blocks have 
dimensions of L x L = (N/M b ) x (N/M b ) with M b ranging 
from 11 to 20. 



